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We consider a family of percolation models in which geometry and connectivity are defined by 
two independent random processes. Such models merge characteristics of discrete and continuous 
percolation. We develop an algorithm allowing effective computation of both universal and model- 
specific percolation quantities in the case when both random processes are Poisson processes. The 
algorithm extends percolation algorithm by Newman and Ziff (M.E.J. Newman and R.M. Ziff, Phys 
Rev E, 64(1):016706, 2001) to handle inhomogeneous lattices. In particular, we use the proposed 
method to compute critical exponents and cluster density distribution in two and three dimensions 
for the model of parallel random tubes connected randomly by bonds, which models the connectivity 
properties of activated carbon. 


I. INTRODUCTION 

Two basic types of percolation models are discrete and 
continuous percolationii^. In the discrete case, a lattice 
is given and its bonds (edges) are open, or its sites (ver¬ 
tices) are occupied, with a probability p, which is the 
relevant parameter of the model. Depending on the case, 
we speak of bond percolation or site percolation. The lo¬ 
cal random variables, which determine bond openness or 
site occupations, define global connections and the main 
focus of the theory is the phenomenon of percolation, i.e. 
the appearance of an infinite cluster (or, in some models: 
of infinite clusters) of connected bonds or sites. 

In continuum models the positions of percolating ob¬ 
jects themselves are chosen at random in space and the 
connections are determined solely by the realization of 
the objects^. A parameter r] playing a role analogous to 
p is usually defined as the expected value of the local den¬ 
sity of the objects. We will usually refer to r] or p as the 


model parameters. In the discrete approach, one can also 
generate the lattice randomly, and then open its edges 
with the same probability, independently of the random 
geometry. Classical examples of discrete and continuum 
percolation are presented in Fig. [T] However, there are 
instances when complexities of percolation phenomena 
are beyond the scope of these two basic types of perco¬ 
lation model. A simple example is a system of roads, in 
which width of a road is describe by the weight of the 
corresponding edge and the traffic intensity corresponds 
to the percolation parameter. In this situation the prob¬ 
ability of a road connection between two points being 
open is a function of both these parameters^. Another 
interesting case, so-called radio tower model^, is obtained 
by modifying the disc percolation mode£. In this model 
we first randomly distribute in the plane points (towers) 
which are the centers of discs with fixed radius R. The 
different towers cannot communicate beyond the distance 
R, which is the parameter of the model. We set the prob- 
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a) discrete percolation b) continuum percolation 

FIG. 1. Examples of percolation discrete and continuum per¬ 
colation: a) bond percolation on the square lattice, b) discs 
in the plane. Clusters are delineated in both cases. 


ability that a connection (an open bond) exists between 
a pair of towers as pbond = niax(0,1 — d/R), where d is 
a distance between two points. We look for the critical 
value of R at which an infinite cluster appears. 

These models have two things in common: their ge¬ 
ometry is random and the possible connections in the 
system are determined by a random variable, whose dis¬ 
tribution is defined by both the geometry and the model 
parameter. Models of such discrete-continuous class are 
thus described by a random graph with weighted edges, 
where the probability of a connection depends both on 
the percolation parameter, and on the weight of edges, 
dictated by the geometry of the graph realization. 

In the present paper we introduce a class of discrete- 
continuous percolation models, consisting of parallel ran¬ 
dom tubes connected randomly by bonds. These models 
merge the characteristics of both discrete and continuous 
types of percolation, and are proposed to describe some 
connectivity properties of activated carbon^. In order 
to compute such properties, however, efficient algorithms 
for discrete-continuous types of percolation models have 
to be developed. To meet this challenge, in particular 
to handle inhomogeneous lattices, we extend an efficient 
percolation algorithm by Newman and ZiffiS. 

Motivations for using such an inhomogeneous tube- 
based model to simulate percolative properties of ac¬ 
tivated carbon originate from wood processing science. 
In the process of wood gasification the material is first 
transformed into charcoal containing approximately one- 
third of its initial mass, and then into various stages 
of activated carbon. Finally the structure of the ma¬ 
terial breaks down making the material collapse into 
fine dust, which burns into a small amount of as h^^d^ . 
Fragmentation, which is a phenomenon closely related 
to percolation, is observed during thermal conversion of 
charcoa l^^d^ . In this process the initial structure (skele¬ 
ton) of wood, composed of parallel cylinders, persists, but 
the hemicellulose, cellulose and lignin that form walls of 


the cylinders are transformed into more carbonic com¬ 
pounds. In this way, although the initial skeleton per¬ 
sists, the microscopic structure of its walls becomes much 
more complex. Fine micro-porous substructures^^ are 
formed, and lead to a rapid increase of internal sur¬ 
face (specihc surface area), when charcoal is transformed 
into activated carbon. Several models were developed 
to explain the complicated micro-structures observed in 
charcoal and activated carbonic. In particular, vari¬ 
ous forms of carbon potentially building such micro¬ 
structures were considered: graphene ribbons (model 
of Jenkins-Kawamura), fullerenes (Harris mode b^d^ ), 
stacked graphiteii or graphenei^, carbon onionsi^ and 
nanotubes^‘2i‘I^. In the present paper we explore perco¬ 
lation properties of a tube-based model, representing a 
nanopipe micro-structure of activated carbon^i^ d^d° . We 
assume that the skeleton walls are made of a collection of 
parallel tubes representing nanopipes of varying lengths. 
These nanopipes form an inhomogeneous lattice bound 
together by amorphous carbon connections. We assume 
that during gasification with CO 2 and H 2 O, the amor¬ 
phous carbon is reacting with these gasification agents, 
and the bonds are removed. The bond removal leaves 
more and more nanopipes disconnected, leading to disin¬ 
tegration of small clusters and, finally, to the breakdown 
on the percolating skeleton. Potential applications of the 
introduced class of percolation models and the developed 
algorithm are beyond this particular tube-based descrip¬ 
tion of activated carbon, including also above mentioned 
models of road networks and radio towers, as well as other 
discrete-continuous percolation systems. 

In Section El we describe a tube-based percolation 
model. We then introduce an extension of the well-knows 
percolation algorithm by Newman and Ziffi^, which al¬ 
lows us to treat the inhomogeneity of the lattice inherent 
in our model fSection lHll) . We validate the extended algo¬ 
rithm in Section llVl comparing its results with the known 
exact solutions for two-dimensional percolation, and use 
it to obtain new results for the three dimensional problem 
in the final Section [V] 


II. PERCOLATION MODEL 

Here we present the tube-based model. To dehne the 
model precisely in two-dimensions, we proceed in three 
steps: 

• we start from n parallel (vertical, for definiteness) 
lines of length L. 

• we use n independent Poisson processes with the 
same parameter /ii to divide the lines into seg¬ 
ments, called tubes. 

• we introduce bonds between each pair of adjacent 
lines and in this manner the connections between 
tubes are established. The bonds are generated by 
independent Poisson processes with parameter ^ 2 - 
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The resulting graph is presented in Fig. [2j 


Breaks, 
generated by 


Bonds, 
generated by 
process 



FIG. 2. An example of a realization of the two dimensional 
graph described by a set of four parameters {L, n, pi, P 2 }. 


A three dimensional model is defined similarly. First, 
we introduce a set of lines of length L passing through 
the points of a square lattice and perpendicular to plane 
of this lattice. Then we follow the procedure for 2D case, 
dividing lines into segments and generating bonds be¬ 
tween each pair of adjacent lines. “Adjacent” is defined 
here using the nearest-neighbor connections on the un¬ 
derlying square lattice, so that in the 2D case a line not 
lying on the boundary has two adjacent lines while in the 
3D case it has four. 

The resulting discrete-continuous model consists of 
parallel tubes of random length connected randomly by 
bonds whose distribution is defined by the spatial loca¬ 
tion of the tubes and by the model parameter. 

The resulting random graph model is described by four 
parameters {L,n, ^ 1 ,^ 2 }; defining the size of the model 
{L and n), length of the tubes described by Poisson pro¬ 
cesses with the parameter pi and with bonds between 
these tubes generated by independent Poisson processes 
with the parameter p 2 - Under rescaling in the direction 
of the lines, the resulting graph is equivalent to the sys¬ 
tem with parameters {Lpi,n, 1,p 2 //ii}- We thus put 
Pi := 1 and p 2 = p, so in the limit when L and n go to 
infinity at the same rate the model has only one param¬ 
eter p. For simplicity in most of the simulations we put 
L = n. 

By definition, different segments (tubes) of the same 
line are not connected to each other. Only tubes lying 
on adjacent lines may be connected, if one or more open 
bonds between them are established. A single open bond 
is sufficient to connect two tubes. This allows one to 
calculate a connection probability between two adjacent 
tubes in terms of their relative position as follows. Two 


tubes lying on adjacent lines may only be connected if 
there is a nonzero overlap h between their vertical po¬ 
sitions as shown in Fig. [31 The probability that two 
such tubes have k open bonds is given by the Poisson 
distribution with parameter fih. That is, 


m 


fc! 


( 1 ) 


Tubes are disconnected (fc = 0) with probability P{0) = 
they are connected with probability 

I^bond ~ 1 6 ^ • 


Possible Overlapping 

connections distances h. 



FIG. 3. Redefined graph. The overlap between two adjacent 
tunes is shown by intervals between the arrows. 

A sample realization of the two dimensional model is 
presented in Fig. [H Groups of connected tubes form 
clusters marked in Fig. 2] by a single color. 


III. THE ALGORITHM 
A. Percolation threshold 

For any percolation model on a square lattice L x L 
one defines the crossing probability n(p, L) as the prob¬ 
ability that there is an open connection between the left 
boundary and the right boundary. The crossing proba¬ 
bility depends on the size of the lattice and on the model 
parameter p. In the limit L —^ 00 , 11 converges to 0 for 
p < Pc and to 1 when p > pc- The critical value Pc is 
called the percolation threshold or the critical point, and 
depends on the type of lattice (e.g. square, triangular, 
etc<^). For a finite lattice, the transition is not sharp and 
many approximations of the critical point are used. Ex¬ 
amples are the point pci, where the crossing probability 
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FIG. 4. A sample realization of tube-based model with clus¬ 
ters of connected tubes marked by a single color: the brighter 
the color the larger the cluster size. 

is equal to the point where the slope of 11 (as a 

function of p) is largest, or, as used in this paper, 

Pav(i) = J p^{p,L)dp. (2) 

The term ^{p,L)dp can be interpreted as the probabil¬ 
ity that the graph begins to percolate for a value of the 
model parameter in the interval (p, p + dp). Thus Pav 
is the expected value of p at the onset of percolation^^. 
Similarly, a measure of the width of the transition region 
can be defined as the variance 

= J (P-Pav)^ ^{p,L)dp (3) 

These quantities satisfy the scaling relations^: 

_ 

Pav - Pc oc L (4a) 

A oc (4b) 

where v is the (universal) critical length exponent. For 
additional discussion see als o^^'^^ . This leads to an 
asymptotic linear relation between p^viL) and A(L): 

Pav = aA+ Pc, (5) 

where a is a proportionality constant. Equation [5] pro¬ 
vides a simple method of extrapolating results obtained 
for finite lattices to the infinite one. 


distribution on the interval [0,1]. To simulate a realiza¬ 
tion with density p, we open the bonds for which < p. 
We then check for existence of an open connection be¬ 
tween the opposite sides of the lattice. Applying this with 
different p (for the same realization of the r^), we approx¬ 
imate Peon as the value of p at which the connection first 
forms for a given realization. The consecutive values of 
p are selected as in the binary search algorithm. Repeat¬ 
ing the whole procedure many times for different sets of 
random numbers r^, we obtain a set of values Pcon- This 
allows us to estimate pav by the empirical mean value of 
Peon and A as its empirical variance. Such procedure is 
the basis of many more advanced methods of computing 
the percolation threshold, such as Hoshen-Kopelman^^ 
and Leath-Alexandrowica^ii^ algorithms. We propose 
to follow a different approach, which is a modification of 
the Newman-Ziff algorithn>iS, computing the value of p 
at which an open connection appears for a given realiza¬ 
tion in a single run. Unlike in the original Newman-Ziff 
approach which used ‘micro-canonical ensemble— we 
use ‘canonical ensemble’. The main advantage of the 
modified approach is its applicability to more general 
graphs, where probabilities vary from bond to bond. We 
note that the transformation between ‘micro-canonical’ 
and ‘canonical’ ensemble representations is complicated 
and impractical in this generality. From the point of view 
of computing percolation threshold on homogeneous lat¬ 
tices both algorithms are equivalent, as explained in de¬ 
tail below. 

The idea of the so-called Rising Water algorithm, in¬ 
spired by a remark iiti is again to assign a random num¬ 
ber to every bond, as described above. To determine 
the value of p at which percolation sets in, consecutive 
bonds i are open in the order of the increasing Vi. As¬ 
suming that random numbers assigned to different bonds 
are different, at each stage we obtain the same graph as 
when using the simplest method described above with 
p = Ti. The algorithm stops when a connection linking 
a fixed pair of opposite sides of the square is established. 
The estimate of Pcon is equal to the value of the last 
added bond. The results of applying the two algorithms 
are identical. Indeed, the first algorithm applied with 
P A Pcon, where Pcon is a result of the Rising Water algo¬ 
rithm, and the same sequence of will find a connection. 
On the other hand, for (p < Pcon) no connection will be 
found, which shows that the two algorithms indeed yield 
the same result. 


C. Extension of the algorithm 


B. Algorithms for the homogeneous lattice 

For the simplest example of an algorithm computing 
the critical density, consider the bond percolation model 
on a regular, homogeneous lattice. We assign to each 
bond i a random number Vi, sampled from the uniform 


In case of the general model studied here, in which 
probabilities of connections depend on both the geometry 
and the model parameter, both Newman-Ziff and Rising 
Water algorithms need further modifications. 

In the simplest algorithm applied to the tube-based 
model we have to generate random values for all pairs 
of adjacent tubes and connect a pair of tubes with an 
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overlap hi when the following condition is fulfilled: 

n > e-^^. (6) 

To define an extension of the Rising Water Algorithm 
we have to compute, for every bond i, the smallest value 
of model parameter mu for which the equation [5] is ful¬ 
filled. We denote this value Hi, thus 

^i =-expri/hi. (7) 

When ^ < fii i-th bond is closed, and when /j, > i-th 
bond is open. 

Then, as in the homogeneous case we sort the set of 
fii in the increasing order and we open the bonds in the 
graph in this order. We estimate the critical value of 
parameter called /Zcon by the first value of p.i at which a 
connection between two fixed opposite sides of a square 
forms. 

The algorithm introduced by Newman and Ziff and the 
extended algorithm proposed in this paper are summa¬ 
rized in Table HI 

The important parts of these algorithms are two main 
operations: 

• finding the cluster containing a given site; 

• connecting two clusters. 

To make these operations efficient Newman and Ziff 
have proposed to represent the connections within a 
graph by a so-called “union-find” (or “disjoint-set”) data 
structure^!. It stores information about connections in 
the form of trees where every site points either to an¬ 
other site from the same cluster, or to itself. The element 
pointing to itself, is the root of the tree and provides the 
cluster’s identihcation. To find the cluster containing a 
given site, we follow the path indicated by the point¬ 
ers until we reach the root. If for two sites we get the 
same root, both sites belong to the same cluster. To con¬ 
nect two different clusters we add a pointer between their 
roots. Two main modifications are commonly used. The 
first one is to always point from the smaller tree to the 
bigger one (“balancing”). It requires storing the infor¬ 
mation about each cluster’s size. The second is called 
path compression: having found the root of an element’s 
cluster, we re-track the path from the element to the root 
again, changing the parent of each site along the way to 
the root. Using such union-find data structure makes 
operations of adding an edge and checking whether two 
sites belong to the same cluster very fast. 

Beside pointer to the parent and size of the subtree, 
one can store additional information in each site’s record, 
such as moment of inertia, position, or the information 
about the cluster’s connection to boundaries. The last 
one is a simple way to check for whether the opposite 
parts of the boundary are connected. 

The position of a site can be used to check whether the 
cluster is wrapped around the torus^i^. 

The amortized computational cost of using it is pro¬ 
portional to the inverse Ackermann function and thus 


it can be considered as a small constant for practical 
purposes^!. 


D. Percolation statistics 

In contrast to older approaches, the important nov¬ 
elty of Newman-Ziff algorith m^^i^^ , as shown in Table U 
(step 4), is its ability to simultaneously calculate a model 
characteristic of a given configuration for different val¬ 
ues of the model parameter p. While standard methods 
need K runs of the algorithm to compute K values of a 
model characteristic for a given set of model parameters 
Pk (fc = 1, ■■■K), in our approach, as in that of Newman 
and Ziff, all values values are obtained simultaneously in 
a single run. Both methods can obtain many important 
characteristics of the model, for example average clus¬ 
ter size, average moment of inertia and so on, with con¬ 
stant computational cost in every run of the algorithm. 
Other parameters like histogram of cluster-size distribu¬ 
tion with B bins can be calculated with an additional 
cost proportional to the number of points in the realiza¬ 
tion (N) and to the number of bins. Let us consider a 
quantity Q. According to the Newman-Ziff algorithm we 
calculate Q[i\ which is a value of Q after adding the i-th 
bond. The values Q[i] are then averaged over K differ¬ 
ent realizations, where the value of K depends on the 
required accuracy. As the next step we transform the 
result to the canonical value Q{p) using Eqn. [5] In the 
Rising Water algorithm we calculate Q[j], the values of 
Q for a chosen collection of values of model parameters 
Pj and take the Q[j] obtained in the last step of the al¬ 
gorithm (described in Table U as a step 4) for which we 
had p < Pi- In our method the possibility of effectively 
achieving statistics is related to the operation on clusters. 
Efficiency of our algorithm relies on fast updates of Q, 
using operations on clusters rather than having to run 
through the whole graph at each step. 

For example, we consider the cluster size. The size of 
cluster C (sc) obtained as a union of two clusters A and 
B is equal to: 


sc = SA + SB (9) 

similarly for the calculation of the moment of inertia for 
clusters we use the stored quantities: sizes of clusters Si, 
masses of clusters mi, centers of mass and previous 
moments of inertia li. For unions of clusters we obtain: 


Wc = ma + mb (lOa) 


raiTia + nmb 


Ic = Ia + Ib + {ra - Tcf *ma + {n “ r^f' * mb 


(10b) 

(10c) 


Note that Eqn. llOcI is the parallel axis theorem (Steiner 
law). In our method, if we store in memory informa¬ 
tion about the clusters, all these operations have only a 
constant cost per operation. For example to get a mean 
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TABLE 1. Comparison of the Newman-Ziff algorithm with an extended algorithm. 


Newman-Ziff 


extended algorithm 


1. create a table C[1 : N] to store statistic 

2. run K times for k=l:K 

a) create a list of all bonds 

b) generate a permutation of connections: ji means that j-th bond 
will be added in i-th step 

c) initialize the list of clusters so that each site is an a cluster of 
exactly one site 

d) for i=l:N do 

- look at bond ji connecting sites a and b. If these sites belong to 
different clusters A and B, merge both clusters 

- check for spanning: for the first occurrence save iteration number 
i as ik 

- refresh the statistics in merged cluster and table Q[i] 

3. compute the percolation threshold using the values oi 


4. compute the transformation from microcanonical Q[n] to 
canonical Q{p) using the following formula 

^ N 

n = 0 


1. for a given set of values of model parameter pf (where I = 

1. 2,...) create a table Q[...] to store statistic 

2. run K times for k=l:K 

a) create a list of all bonds 

b) assign a random number ri to every connection and compute 
value of model parameter pi {pi from Eqn. [7|l for which we add 
the bond. Sort connections in order of increasing pi. Let ji denote 
a sorting permutation 

c) initialize the list of clusters so that each site is an a cluster of 
exactly one site 

d) for i=l:N do 

- look at bond ji connecting sites a and b. If these sites belong to 
different clusters A and B, merge both clusters 

- check for spanning, for the first occurrence save pi number as 
Pcon,k 

- refresh the statistics in merged cluster and if for any i, Pi—i < 
pf < Pi, update the statistics Q\pf] 

3. compute the percolation threshold pav and its variance Aav 

using Pcon,k as follows: Pav = = l Pcon,k and Aav = 

j ^k — liPcon,k ~ Pav)^ 


value of the moment of inertia we additionally store in 
memory the sum of the moments of inertia of the clusters 
and update this sum. 


E. Critical exponents 

When the percolation threshold pc is computed, a post¬ 
processing algorithm gathers statistics about the distri¬ 
bution of clusters (including the size of the largest cluster, 
cluster-size moments, cluster-volume moments). These 
statistics are determined for p in a vicinity of Pc- This 
allows computing several critical exponents of the model. 
In particular the cluster-size distribution near the perco¬ 
lation threshold allows to compute the Fisher exponent r. 
The P exponent is computed from the size of the maximal 
cluster. From data acquired in the algorithm outlined in 
Sec. IIII Cl exponent u in Eqn. I4bh an be computed using 
the scaling relation (Eqn. I4bl) . 


IV. RESULTS IN THE TWO-DIMENSIONAL 
CASE 

A. Percolation threshold 

The simulation was run for several square lattices with 
size ranging from L = 200 to L = 10000. The estimators 
of Pav and A were acquired for mutually perpendicular 


directions, denoted by NS (top to bottom) and WE (left 
to right). The percolation threshold for the infinite lat¬ 
tice {L —>■ oo) was computed by fitting the data to the 
scaling properties described by Eqn. [5] as presented in 
Fig. [S] The results for the infinite lattice based on the 
intercept of the fitted linear function are the following: 

Pc NS = 0.99999 ±2.5 X 10“® (11a) 

Pc WE = 0.99999 ± 5.0 X 10"® (lib) 

It is worth noting that values Pav converge to Pc from 
both directions, as presented in Fig. [5l The obtained 
value of Pc equal 1 is clearly model-specific, as discussed 
in Section IIVBI 


B. Duality and exact analytic result 

We consider a realization of the two dimensional graph 
defined by {L,n,pi,p 2 } presented in Fig. [6^. We define 
the graph dual to the initial one according to the follow¬ 
ing procedure: 

• dual lines are introduced, each line is placed be¬ 
tween two existing lines; 

• dual lines are divided into tubes (dual tubes) by the 
bonds of initial graph (vertical segments marked in 

Fig. [5b); 
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2D case 
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New tubes and bonds are generated in the same way 
as the original ones, with the two Poisson process param¬ 
eters interchanged. Notice that the two graphs have no 
intersections. We either have a connection from top to 
bottom, using tubes and bonds of the original graph, or 
we can draw a line through the empty spaces and breaks 
between the tubes from left to right, that does not cross 
any bonds or tubes. In the latter case, there is a con¬ 
nection from left to right in the dual graph. Similarly, 
exactly one of the two alternatives occurs: either there is 
a connection from left to right by bonds and tubes of the 
original graph, or there is a connection from top to bot¬ 
tom in the dual graph— an unbroken path through empty 
spaces. A given realization starts to percolate when the 
dual graph stops percolating, so Pav = Pav'^* fo'' ^ of 
dual graphs. 

We know that the percolation threshold in the limit 
n = A —> oo depends only on the ratio P2 /ai- Increas¬ 
ing /ii results in more (shorter) tubes and thus makes 
percolation more difficult, while increasing /i 2 makes for 
more connections between tubes, which facilitates it. To¬ 
gether with the duality described above, this indicates 
that ^ = 1, i.e. pc = 1 is the percolation threshold, 
thus explaining the numerical result (lllal) and dill) , and 
giving further support to our method. We emphasize 
that a rigorous proof that the critical value of p equals 1 
requires a more careful argument. The first result of this 
type (for the square lattice) was proven in2^. Simpler 
arguments developed later can be found in^. They can 
be adapted to cover the present case as well. 


FIG. 5. The percolation threshold computed by studying top- 
to-bottom (green points) and left-to-right (red points) connec¬ 
tions. In b) the differences between data points and the fitted „ „ . . , 

, C. Critical exponents 

line are shown. ^ 


• at the positions on breaks between initial tubes the 
dual bonds connecting dual tubes are introduced 
(horizontal lines marked in Fig. |BJd). 

The two graphs, initial and dual, are shown in Fig [5^ 
and c. 



FIG. 6. Gonstruction of dual graph from the original one. 


Based on the scaling law (Eqn. 0^ and I4b|) we obtain 
the correlation length exponent : ly = 1.345 ±0.009. The 
exact value is known to be 4/3. 

We determined two characteristics of the clusters: the 
first one, presented in Fig. [7]a), based on size of clusters 
and the second one, presented in Fig. [7]b), based on 
volume of clusters. 

The Fisher exponent r, is determined based on clus¬ 
ter size distribution presented in Fig. [7] a) as r = 
2.046 ± 0.023. The exact value is 187/91 « 2.054^. The 
agreement of the results with the known values of critical 
exponents supports the validity of the algorithm. 

Moreover, we show that the slopes of lines fitted in 
Figs. [7] a) and b) are the same, thus the Fisher exponent 
determined based on cluster size distribution and the ex¬ 
ponent which based on cluster volume distribution are 
also the same. This observation confirms the duality re¬ 
lation of percolation models on a given and dual graphs, 
discussed in Section HVB I 




























































































b) Normalized number of clusters of volume v. 


FIG. 7. Number of clusters: a) ris - number of clusters of size 
s per one site, b) Uv - number of clusters of volume v per unit 
volume. Data from 2D grid 15000 x 15000. 


V. RESULTS IN THREE DIMENSIONS 

The simulation was run for cubic lattices with size 
ranging from L = 100 to L = 400. The estimators of 
Pav and A were acquired for perpendicular directions, 
denoted by NS, WE and TB (top to bottom). As in 
m we use scaling properties described by the Eqn. 133 
to compute the percolation threshold for infinite lattice. 
The results are as follows: 

p^NS = 0.231466±6 X 10“® (12a) 

Pcwe = 0.23146±7x 10“® (12b) 

Pc TB = 0.23140 ± 1.2 X lO"'^ (12c) 

The results obtained by fitting independently three linear 
functions, as presented in[8]a), can be improved using the 
following constraints: 


• the line fitted to the results parallel to tubes (TB) 
have the same intercept b. 

Thus the improved estimated value of the percolation 
threshold is: 


Pc = 0.231456 ± 6 x 10“® (13a) 

It is worth noting that, exactly as in the two-dimensional 
case, which we discussed in Section (lIVBjl . the values pav 
converge to pc from both directions. This is clearly visible 
in Fig. [H] 


3D case 



a) 


15 
10 
5 

of 0 

-5 
-10 
-15 

0 0.001 0.002 0.003 0.004 0.005 

A 

b) 


3D case 



FIG. 8. The percolation threshold computed from top to bot¬ 
tom connections (green points) and from left to right connec¬ 
tions (red points). In b) the differences between data points 
and the htted lines are shown. 


VI. DISCUSSION 


A. Computational cost 


• the lines fitted to the results perpendicular to tubes 
(NS and WE) have the same slope and intercept b; 


The computational cost of determining an approximate 
value of the critical point pav depends on the size of the 
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lattice and on the desired accuracy of calculation which 
can be expressed in terms of standard deviation A(L). 
Analysis of this computational complexity allows us to 
know what accuracy e can be achieved in a given time. 
The obtained value of Pav is approximated by the Monte 
Carlo estimator pav, which takes into account all runs of 
the algorithm: 

Apav CX ^ (14) 

where k denotes the number of repetitions of the Rising 
Water algorithm. Thus e, the final accuracy of Pav, de¬ 
pends on the number of runs of the algorithm and on the 
variance A(L) as follows: 


A(T) 

\/fc 


(15) 


From the Eqn. |4b]we know that A(L) depends on the size 
of the domain L. The computational cost c = kL‘^ log n 
is proportional to the number k of times the Rising Water 
algorithm is repeated and to the cost of a single run (of 
the order of L‘^ log n, where d is the dimensionality of 
the problem). Thus the computational cost to obtain 
the result with the accuracy e is 


c = 



(16) 


The exponent 2 — depends on the dimension of the 
problem. For the two-dimensional case it is 1/2, while in 
three dimensions it equals approximately 0.72. The log¬ 
arithmic factor in the expression for the computational 
cost (Eqn. I16|) is due to sorting of random numbers in 
step 2c in Section HI One method to avoid this is to use 
so-called bucket sort, which is a linear-time sorting algo¬ 
rithm using information about data distributional. Due 
to statistical behavior of the random values pi we can 
create a set of disjoint intervals that cover all possible 
values of pi and have approximately the same expected 
number of random values pi in each interval. Let us de¬ 
note this expected number of random variables in one 
interval (“bucket”) by M. For every generated random 
Pi {i = 1,...,A^) we can compute in constant time to 
which bucket it should be assigned. When all numbers 
are generated and classified, in each bucket we have a 
set oi M + O numbers, and we need to sort it. 

The computational cost of generating N random vari¬ 
ables and sorting N/M buckets of size M is 0{N log M) 
and it is linear in N because it is always possible to gen¬ 
erate enough intervals to keep M constant. After that, 
the cost of running the algorithm k times is 

0{kL‘^a{L'^)) 


and cost of running the algorithm to get desired accuracy 
e is 


O 




Here a{L‘^) is the Ackerman function and can be consid¬ 
ered constant. Despite better asymptotic behavior of the 
bucket sort, it does not give a better performance except 
for very big lattices. 


VII. CONCLUSIONS 


In summary, three goals have been achieved in this 
work: 

• We have defined a family of discrete-continuous 
percolation models motivated by the physics 
of activated carbon. These models deal with 
tubes of random length connected by random 
bonds; as such they should describe well situa¬ 
tions in which complicated micro-structures ob¬ 
served in activated carbon have approximately lin¬ 
ear textures: graphene ribbons (model of Jenkins- 
Kawamura) and nanotubesi^. In cases the struc¬ 
tures are neither ID nor quasi-lD (fullerenes (Har¬ 
ris mode b^d^ ), stacked graphiteii, or graphenei^, 
carbon onionsi^), the concrete models considered 
here provide only a “caricature” of the real situa¬ 
tion. Still we expect that even in these cases they 
capture some qualitative aspects of the underlying 
physics. 

• We have extended the standard algorithm of New¬ 
man and Zifii^ to handle inhomogeneous lattices. 
This extension is non-trivial, and we have analyzed 
in detail its convergence properties. 

• We applied the extended algorithm to the family of 
models in question, calculating critical parameters 
and cluster density distributions in two and three 
dimensions. 

Possibilities for further studies include: i) applications 
of the present models to experimental data, suggest¬ 
ing geometry formed by parallel random tubes/ribbons 
connected randomly by bonds; ii) development of con¬ 
crete models with geometry formed by parallel random 
flakes/patches connected randomly by bonds; iii) appli¬ 
cation of the method to such models, calculation of their 
properties, and direct comparison with experiments. 

It is worth mentioning that the problem of quantum 
aspects of the carbon activation process is also to a great 
extent open. This suggests to study quantum versions 
of the family of the discrete-continuous models discussed 
in this paper. The interplay of discrete and continuous 
aspects may lead to quantitatively novel effects. It is 
worth noting that such quantum disordered models can 
in principle be simulated, quantum simulated, by a sys¬ 
tem of ultracold atoms (see, for instance^): an array of 
random length ID Bose condensed gases with controlled 
random connections between them. 
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